{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Identifying the diffusion equation from a random walk\n",
    "\n",
    "Samuel Rudy, 2016\n",
    "\n",
    "Here we take various lengths of a random walk where $x_{j+1} \\sim \\mathcal{N}(x_j, dt)$ and see if we can identify the diffusion equation.  As expected, it works better for longer series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (12,6)\n",
    "import numpy as np\n",
    "from PDE_FIND import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trial: 1 of 10\n",
      "Length of time series used:  100\n",
      "u_t = (-121283487.459491 -0.000000i)u^3\n",
      "   \n",
      "Length of time series used:  146\n",
      "u_t = (-105893062.687726 -0.000000i)u^3\n",
      "   \n",
      "Length of time series used:  215\n",
      "u_t = (-40949708.875765 -0.000000i)u^3\n",
      "   \n",
      "Length of time series used:  316\n",
      "u_t = (-25929330.762661 -0.000000i)u^3\n",
      "   \n",
      "Length of time series used:  464\n",
      "u_t = (-9727002.518794 -0.000000i)u^3\n",
      "   \n",
      "Length of time series used:  681\n",
      "u_t = (0.097842 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  1000\n",
      "u_t = (0.132088 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  1467\n",
      "u_t = (0.292211 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  2154\n",
      "u_t = (0.336024 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  3162\n",
      "u_t = (0.401712 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  4641\n",
      "u_t = (0.409962 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  6812\n",
      "u_t = (0.455302 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  10000\n",
      "u_t = (0.436203 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  14677\n",
      "u_t = (0.461215 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  21544\n",
      "u_t = (0.455262 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  31622\n",
      "u_t = (0.463378 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  46415\n",
      "u_t = (0.470980 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  68129\n",
      "u_t = (0.470179 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  100000\n",
      "u_t = (0.482274 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  146779\n",
      "u_t = (0.475499 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  215443\n",
      "u_t = (0.470402 +0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  316227\n",
      "u_t = (0.475082 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  464158\n",
      "u_t = (0.483958 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  681292\n",
      "u_t = (0.493118 -0.000000i)u_{xx}\n",
      "   \n",
      "Length of time series used:  1000000\n",
      "u_t = (0.493998 -0.000000i)u_{xx}\n",
      "   \n",
      "Trial: 2 of 10\n",
      "Trial: 3 of 10\n",
      "Trial: 4 of 10\n",
      "Trial: 5 of 10\n",
      "Trial: 6 of 10\n",
      "Trial: 7 of 10\n",
      "Trial: 8 of 10\n",
      "Trial: 9 of 10\n",
      "Trial: 10 of 10\n"
     ]
    }
   ],
   "source": [
    "# set seed = 0\n",
    "numpy.random.seed(0)\n",
    "\n",
    "N = 24\n",
    "lengths = [int(10**j) for j in 2+np.arange(N+1)*4.0/N]\n",
    "err = {}\n",
    "for l in range(N+1): err[l] = []\n",
    "sparsity_err = np.zeros(len(lengths))\n",
    "trials = 10\n",
    "w_true = np.zeros((16,1))\n",
    "w_true[8] = 0.5\n",
    "\n",
    "for trial in range(trials):\n",
    "    \n",
    "    print \"Trial:\", trial+1, \"of\", trials\n",
    "    \n",
    "    # generate a new time series\n",
    "    dt = 0.01\n",
    "    advection = 0 # it has some trouble with advection\n",
    "    pos = np.cumsum(np.sqrt(dt)*np.random.randn(lengths[-1])) + dt*advection*np.arange(lengths[-1])\n",
    "\n",
    "    # fit various lengths of it to a pde\n",
    "    for l in range(len(lengths)):\n",
    "\n",
    "        length = lengths[l]\n",
    "        P = {}\n",
    "        M = 0\n",
    "        m = 5\n",
    "        \n",
    "        # More bins for longer time series.  We just need to make sure there aren't too many bins for how many points we have\n",
    "        n = int(20*log(length))\n",
    "        for j in range(m): P[j] = []\n",
    "\n",
    "        for i in range(length-m):\n",
    "\n",
    "            # center\n",
    "            y = pos[i+1:i+m+1] - pos[i]\n",
    "            M = max([M, max(abs(y))])\n",
    "\n",
    "            # add to distribution\n",
    "            for j in range(m):\n",
    "                P[j].append(y[j])\n",
    "\n",
    "        bins = np.linspace(-M,M,n+1)\n",
    "        x = linspace(M*(1/n-1),M*(1-1/n),n)\n",
    "        dx = x[2]-x[1]\n",
    "        \n",
    "        U = np.zeros((n,m))\n",
    "        for i in range(m):\n",
    "            U[:,i] = np.histogram(P[i],bins)[0]/float(dx*(len(pos)-m))\n",
    "\n",
    "        Ut,R,rhs_des = build_linear_system(U, dt, dx, D=3, P=3, time_diff = 'FD', deg_x = 5, width_x = np.max([10,n/10]))\n",
    "        lam = 10**-3*length\n",
    "        d_tol_2 = 0.001/dx\n",
    "        d_tol_0 = 0.001/dx\n",
    "        \n",
    "        # Try two different normalizations and see which one performs better.  They should get the same answer for most of \n",
    "        # the longer runs.\n",
    "        split = 0.8\n",
    "        N = len(Ut)\n",
    "        train = np.random.choice(N, int(N*split), replace = False)\n",
    "        test = [i for i in np.arange(N) if i not in train]\n",
    "        \n",
    "        w2 = TrainSTRidge(R[train,:], Ut[train], lam, d_tol_2, normalize = 2)\n",
    "        w0 = TrainSTRidge(R[train,:], Ut[train], lam, d_tol_0, normalize = 0)\n",
    "                \n",
    "        err2 = np.linalg.norm(Ut[test] - R[test,:].dot(w2)) + 0.01*np.linalg.norm(Ut[test], 2)*np.count_nonzero(w2)\n",
    "        err0 = np.linalg.norm(Ut[test] - R[test,:].dot(w0)) + 0.01*np.linalg.norm(Ut[test], 2)*np.count_nonzero(w0)\n",
    "        \n",
    "        w = [w0,w2]\n",
    "        error = [err0,err2]\n",
    "        method = argmin(error)\n",
    "        w_r = w[method]\n",
    "        err[l].append(np.linalg.norm(w_r - w_true, 1))\n",
    "\n",
    "        if trial == 0:\n",
    "            print \"Length of time series used: \", length\n",
    "#             print \"Condition Number: \", cond(R)\n",
    "#             print \"Regularization: \", ['unregularized','2 norm'][method] \n",
    "            print_pde(w_r, rhs_des)\n",
    "        \n",
    "        w_r = np.array(w_r)\n",
    "        sparsity_err[l] += float(len(np.where(w_r[0:5] != 0)[0]) + len(np.where(w_r[7:8] != 0)[0]) + int(w_r[6] == 0))/trials\n",
    "        \n",
    "# print err\n",
    "# print sparsity_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fb2f790c190>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuQAAAGaCAYAAABHZFZJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X2cXGV58PHfRRAw8hqKChFBQZRCi0SToI2wICj19ZFW\nLUKM1epDC5S39rEVkBAsjy0VSqNQ+whWY3ypyJvaYlFchCiIotSCkQLZFAFBEmTBGAjkev44Z8Nk\nMrvZPTs7Z3fm9/185rOz5z5zn2vOnJ295p77XCcyE0mSJEn12KLuACRJkqReZkIuSZIk1ciEXJIk\nSaqRCbkkSZJUIxNySZIkqUYm5JIkSVKNTMglqYtFxOsi4rsR8UhErI+Iy+uOqVFELCzjOrhp+fqI\nuK7F+s+LiM9ExL0R8VREPB0R25dt20XEP0bEiohYV7b97gTHv0cZ66UTuR1tyn2vbrJl3QFIU1lE\nrG9atB54BPhP4FOZ+YXORzW5lfusPzMPqzuWbhcRewBXUhyTlwCDwPJag9pUlrfRLv8McDjwBeCu\ncp21Zdt5wAeArwKfBZ4GftHmeFsZLlaN0yjeL9z36gom5NL4JbAQCOBZwMuAtwKHRsQrMvMvaoxN\nve1wYGvg1Mz8Ut3BjNG+wJrGBRHxLIrndG1mzm/xmDcCP8vMt3YgviH3UcT6aAe3qYL7Xl3DhFxq\ng8w8p/H3iDgU+CZwckT8Y2b+Tz2RqcfNLH8+UGsUFWTmnS0W70ox1XK457MbcP2EBdVCZj4FtIpV\nE8x9r27iHHJpAmTmtymmBgQwe2h5RLwnIi6LiLsjYk1EPBoRN0bEMa36iYj+ch7ssyLiwxGxPCLW\nDs2ZjIjtI+IvI+Jb5ZzaJyLioYi4KiIOGqbP9RFxXUQ8NyIujYhfRMTjEbEsIuaV60yPiPMiYqDc\n3n9FxB8O93wj4uiI+HY5T/k3EXFHRJweEVs1rLOg/Po5gb4yjqHbh5v6m1vupwfK5/Q/EfFPEbHr\nWPfR5kTESyPiX8ptPFHuj6URsU+Ldf+ljHfPiDgxIm4rX8fryvZDhp5PRMyOiK9HxKoyvhc29DMr\nIr4SEQ+WsQ5ExCci4vlj3eYwz+mQcl8vpDgG+8s+no6GudoRsXdEfDYifl4+9/uimJ+9d4s+N8z1\njoh3RcRNEfFYRNwzyv38ioi4JiIGy+P+2uGO0XL9jeaQR8QKYIDi+HlPw7Hz6fLYG5o+1nhsDb0u\nQ+u/ezTbKpdtGxFnRsRPyngHI+KuiPhiRBzYsN6w85gj4vnl67oinvnb/EpEzGqx7oKhGCPi0PI5\nDe2rr0XEy0bcwa2f1+sj4t8i4pflcXZXRPxdROwwzPqHR8QNUbwfrIqIKxr+PtY3HcMbjvVh+hpo\nPjZiDO9XMYr3i8m876WxcoRcmjhR/myc33gR8F8Uo3gPADsDbwCWRMQ+mXlWUx9Dj/0K8Erg34Er\ngIfK5fsCHyn7+xrFXOEXAm8Bfj8i3pSZ/9Eith2BZRRzij8PzACOBq6JiFcD/1yu81WKaThHA1+M\niFdn5vc3epLFP8P3APcClwG/Ag4CzgEOi4gjMnM98COKBHEhRWL1Lw3d9Df0917gkxTzgq8u+30J\n8D7gzRExNzN/PoZ9NKyIOLJ83Jblc70LeAFwFPDGiOjLzB83bSuBfwTmAV8vb083df1q4EPADRRz\nt38LeLLc5pvK/UT5cyXwCuBPgbdExLzMXFlhm40GKPZzH3AIxbzrgYY2ImI2xbc4z6HYz3dQTLc6\nFnhrRLw2M3/YIo6/oJg28lXgOqBlcteoPKaupTiWvgLcDbyc4nUf9oNFkwuAPYGTgR9TzI2nvL8j\n8G02PbYGmmIfi28ArwK+C/w/4CmKY+NQ4DsUx/OwImJPir+x51M8x88DuwNvpzi2jsrMf2t6WAJv\nppjy9m/AxcBvU0zFeWVE/HZmrh5N8BFxFnAWsIriveEh4HcpXr/fj4hXZebjDev/IfBF4Iny5y8o\njrfvUZwTM9b912r9sbxfjer9opW6971USWZ68+at4o3iJM6nWyw/nCJhegrYvWH5i1qsuyVFYvQE\nsGtT27fLbfwY2KnFY7cDZrRYvhvF/Mrbh4sZ+ETT8mPLtlUUyc5WDW3zyravND3mPeXyLzeuX7Z9\nuNzOiS22f90w+/Ml5X74GfD8prZDy/3ZHMOI+2iE125HioTgQeClTW2/DTwG/KBp+afLbd0LvLBF\nn4c07N8/adH+nHL/rgNe3dT2l+VjrxnLNjfzHM8qYzm4RdtPy7Y/alr+9nJ7d7Toa325X353jHEs\nL7f1pqblJzbsr4Ob2jY5ToA9yuWXjvD3uMmxBSwot/Hu0TwO2L9cdtkw6++wuZgoEvqngb9qWn5Q\n+fr/EpjeFON6ig9ufU2PObfs6y9Gub8PLfu6Adiuqe3dZdvHWhyXTwAHNq3/sYbX6IUNy4eO9Q8P\nE8MK4J6mZVXfr4Z7v5h0+96bt6o3p6xIbRARZ5W3j0TEZRSjtAAXZOa9Q+tl5ormx2YxD/ITFIn5\na1t0n8AZmflIi8c+li1GbTLzforR15dFxAta9LkG+D9Nyz5PkfDuCJyUmU829HcjxSjVy5secxLF\nP7j3Na5f+giwGmg5HWcYf0axH07OzI2qY2QxDehqilHy5zQ9bth9NIIFwPbAwsz8WdO27qAYFT2w\nxdfVCfxtjnxewI8z81Mtlr8V2An4YmZ+t6ntfIp9fESL12w02xy1csT6pcB3M/OLG20o88vAjcBL\no5zC1OSTmfmfY9zWPsD1mfm1puZPUIyWT1ZrWy3MzBFPIoyImcARwP9QVH5pfOxNFBViZlB8E9Ps\nC5nZ37Tsnym+cZszqqjhzymOmQ9k5mNN2/8sxYfXxr/LoeNyaWY2j/yfTZtOmhzH+9WoTYJ9L1Xi\nlBWpPYbmUSbFlI3rgUuyqexhROwO/BVwGMVXtc9uaE6eOQmv2S3DbTgifo8iMT4IeC6wVUPzUJ8/\nb3rYnZn568YFmbk+Ih6kGDlayabuo+GfUkQ8m+Ir8F8Cp0RE8/pBMeK273CxtzA0j7QvIlr9A3wu\nMI0iwWtOHIbdR5vZ1svLr/ebDc0h35dNSwVublvfH2b5LIrX5NvNDZn5dER8B5gPHMimr9lYn99I\nhubRbhJH6Trg98o4bmwMs0IcQ9v6TnNDeczdCLx4jH1OtDsoktajy+kPV1Hshx9k5rpRPH5ojvkN\nmdlqatF1FN9IHQh8rqnth5uuztCH+p1GsW14ZiT4HS3+LqF4j9glInYqP8QOHZetXqPBiPgxcHBz\nWxUV36/Gou59L1ViQi61QWZO29w6EfEiimRmB4qvkr9BMfL0NMXc2AUUJepa9f/gMH2+jWK6yG8o\n5ujeDfya4uvXQyn+ibbqc7gRr6c209b4nrETRdK9C898IGkZ/ghtzXYuf45UKjKBbTdZOMw+2sy2\nAviTzay3ybbYfG3r4dqH5lsPVyVkaPmOFbY5FjtQ7MeR4og2xTG0reFen07UCR+T8oPCoRTH9R8C\nH6XYH49FxGeAv27+QNuk6us89IG+OZ6ny8R6s+8zpZ3LdTf3d7ktxbStoXgn9DUax/vVWNS976VK\nTMilzjmNIol9T2YuaWyIiD+imI89VudQjEK/IpvKxEXEbrRpVGsYQ4n7jzLzlW3uc/vNJDzt2lZS\nzIe+fYyP3dyHjOHah57fJtVUSrs2rTeWbY7FoxQJ5khxZJviGNrW84ZpHy6GdlpfxrDJ/7zhKo6U\n01JOA06LiBdTzJn+38AJFEnfghG2N57XuR0eBSIzf2sM68PYXqOhqjbD5RFD52g06sT7Vd37XqrE\nOeRS5+xV/mx16fI+qiVce1GcfNf8zy2A11Tob9TKhPl2YL+IaDWSOpz1DD/adFP5cyI/SDRuKzq0\nrSE/KrfZ19wQEdN45jW7tQNx0CqO0tBVEdsRx1AfhzQ3RMQWFCcMT7ShxHD3Fm2zWyzbSGbek5mf\npthfj1PMuR7J0P6dVz7HZodR/L1P1Ot8E7BTRIx2utitFMdlq9doezY9dwRG2KdRlM1s9UGnyvvV\nSO8XrdS976VKTMilzhkof/Y1LoyI11OU9Kva50ti0/rVZzO2udtVnU/xFfOnW400RsSO0VCzubSK\n1okRwMcppsZcEBEvadHfs4Y50bCKT1N8RX1WWQKweVsREZskKON0JcWJrkdHxNymtlOAF1FchXI8\nc2g3KzOXUVSymRcRf9DYVpa/m0dxxcsbWz1+jNv6brmtgyPiLU3NJ/LMB9WJ9AOKxO5d5bkPAETE\nDOBvafowHEXN9xe16GcGxfG+pkXbBpl5H8WUjD0pXtfGvudSlBFdTVGecyJcQJFg/79oXbt/etPx\ndxVFgv2uiHhF0+pn0zq5Xk5RNvWtEbFhJD4itqEo0dnKAGN/vxrp/WITk2DfS5U4ZUXqnIuAPwYu\nKyux3E9RXu31wL8Cf1Shzwso6uX+OCK+QnEi1+9R/HO7mqKu7oTJzE+XF9r4M+DuiPgGRXWDGRTJ\n5cHApWX7kG8B74yIqylGqdYB38nMGzLzZ2Ud8kuA2yPiGoor8T2L4iTY11DUU/7tNsS+ukw+Lwdu\niohvUYz4J0UC8KryeUwf77Yatvnr8vn9K3B9RHyZYn+9AngdxTFxXLu2txkLgP8AvhQRV1EkWC+j\nGP19lKI8Xru8r9zWVyLicop67y+nGK38d+DINm5rE5n5i4hYSnEy348j4usUFXbeQHECdvOHxgOA\nyyPiForykPdTnCvxVor/m387is0eR3Ei6N9FxOsoPhS8kGJO+tPAH7eYltXyDMyxyszrIuKDwP8F\n/jsi/o2iDOG2FKUCD6E4j+UN5fq/jogPUNQfvyEivkQx13oesB/FPjq4aRtPRcSFwBkU+/QKin1z\nBMUJ4Pe3CK3K+9Ww7xcj7ILa9r1UlQm5NH6jmmqSmT+JiD6KcoBvoPj7uw14G8VI0zuH6WvY/jPz\nnyNiLcXFUt5NcbLUdyjmo/8hrf/B5WZiHlNbZp4YEf9O8U/wtRRzR1dTJJp/CyxteshJFKOVrwV+\nn+KburMpEgQyc2lZ1eE0ihO9jqA48et+ihPCvjTGmId/MkXiMnSxlNdTJCBPltv6Fs9cwGcs2xpx\n/2bm1WWliQ9RJOE7UJw0dxHwkeZyj6Pc5phl5vfLbwbOoKib/ybgYYrX6yOZ+d9t3NZ3I+I1wN/w\nTPJ9E8W3RUfSOiEfbj+OtH9HavsTiv18NMUHxP8B/oGizvY7mh73A4pk9hCK42InimpCtwD/mJte\nbGuT7Wbmioh4JcX+fUPZ1yDFRWfOzY0vutTYz3A293e78cqZ50XEMooSiPMoLr7zKEWy/E8U5f8a\n1/9KFBfKOouiFv0TFIn4q4C/psXUrsw8KyJ+Dby/vP2i7Pdsig8yzfukyvvViO8XTMJ9L1URmR5j\nkiSptYj4NEUC/aJ21cKXtDHnkEuSJEk1MiGXJEmSamRCLkmSNsf5rdIEcg65JEmSVKOerrISEX4a\nkSRJUkdkZssSmz2dkAMsX7687hDaYvHixZx44olTfpvj7XM0jx8cHGT27I2vA7Nw4UIWLlw4qm2M\nZt12rDOWmCazup5Hu7fbjv6q9NHuY3M063XLsbc53XJstqPPiT42R7v+eNfppmO3jufSLcfmWB/X\nqffO4qK0rU3rlgO3irPPPnvhCSecUHcYbfOCF7ygK7Y53j439/gnnniCmTNnbrJ8zz33HPU2RrPu\neNfp7++nr69v1DFNZmPZt5N5u+3or0of7T42R7NeXa9Zp3XLsdmOPif62Bzt+uNZp5veN6Ge47Nb\njs2xPq4T751nn302CxcuPLtVW0/PIY+I7JYRco1eqxHyyaibRnokqRN839RkFhHDTlmxyoo0SXXT\nKI8kdYLvm5qqHCF3hLznTJURckmS1D0cIZckSZImKRNySZIkqUYm5JIkSVKNTMglSZKkGpmQS5Ik\nSTUyIZckSZJqZEIuSZIk1ciEXJIkSaqRCbkkSZJUIxNySZIkqUY9n5BPW7WKHZcurTsMSZIk9agt\n6w6gbrsvWMDWd90FwK+OOabmaCRJktRren6EfOu77uKJvffmsSOPrDsUSZIk9aCeT8gB7j//fJ7e\neee6w5AkSVIPMiEHdjv1VKatWlV3GJIkSepBPZ+QP7H33mx9111sd801dYciSZKkHtTzJ3Xe+5nP\nsN0113hCpyRJkmrR8yPkT++8s8m4JEmSatPzI+SLFy9mzpw5zJ07t+5QJEmS1GX6+/vp7+8fcZ3I\nzM5EMwlFRC5fvrzuMNRhg4ODzJ49u+4wJElSD4kIMjNatfX8lBVJkiSpTibkkiRJUo1MyCVJkqQa\nmZBLkiRJNTIhlyRJkmpkQi5JkiTVyIRckiRJqpEJuSRJklQjE3JJkiSpRibkkiRJUo1MyCVJkqQa\nmZBLkiRJNTIhlyRJkmpkQi5JkiTVaMu6A5gIEbEF8BlgN2AF8IHMXF9vVJIkSdKmunWE/G3APZn5\nWmA5cFTN8UiSJEktdWtCvhfw4/L+j4CDa4xFkiRJGtakTsgj4viIuCUi1kbEpU1tO0XEFRHxeESs\niIijG5rvAA4r7x8O7NSpmCVJkqSxmNQJOXAfcA5wSYu2i4C1wC7AscDFEbEvQGZ+DVgbEd8EpgO/\n6Ey4kiRJ0thM6oQ8M6/MzKuB1Y3LI2I6xbzwMzLzN5m5DLgKmN/w2L/MzMPLx17VwbAlSZKkUZvU\nCfkI9gHWZebdDctuA/YDiIjnRcR1EXEt8ERm3lhHkJIkSdLmTNWyh9sCg03LBoHtADLzQZ6ZQz6i\n+fPnM3PmTGbOnMmcOXOYO3dueyOVJElSz+nv76e/v5+BgQEGBgZGXHeqJuSPA9s3LdsBeGysHS1Z\nsqQtAUmSJElD+vr66Ovr2/B7RAy77lSdsnInsGVE7NWw7ADg9prikSRJkiqZ1Al5REyLiG2AaRQJ\n+NYRMS0z1wCXA4siYnpEzAPeDDjcLUmSpCllUifkwBnAGuCDwDHl/dPLtuMpSho+BHwOOC4zf1pH\nkJIkSVJVkZl1x1CbiMjly5fXHYY6bHBwkNmzZ9cdhiRJ6iERQWa2nEg+2UfIJUmSpK5mQi5JkiTV\nyIRckiRJqpEJuSRJklQjE3JJkiSpRmNOyCNiq4i4IiIOnoiAJEmSpF4y5oQ8M58EDq/yWEmSJEkb\nq5pULwMOamcgkiRJUi/asuLjTgOujIjHgSuBB4CNrjCUmevHGZskSZLU9aqOkP8E2Au4EFgJPAms\na7g92ZboJEmSpC5XdYR8EU0j4pIkSZLGrlJCnpkL2xyHJEmS1JOslCJJkiTVqHJCHhG7RsTfR8Qt\nEXF3+fPvIuL57QxQkiRJ6maVEvKI2Af4MfDnwOPA98ufJwE/joiXtC1CSZIkqYtVPanzb4FBYG5m\nDgwtjIg9gP8o248ad3SSJElSl6s6ZeVQ4MzGZBwgM1cCC8t2SZIkSZtRNSHfCnhsmLbHynZJkiRJ\nm1E1If8xcGJEbPT4iAjgz8r2KWHx4sXcfPPNdYchSZKkLtTf38/ChQtHXCcyx359n4g4EvgacDfw\nJeAB4PnA24GXAG/MzP8Yc8cdFhG5fPnyusNQhw0ODjJ79uy6w5AkST0kIsjMaNVW9cJA10TEm4CP\nAKcDQXHlzh8Cb5oKybgkSZI0GVSuQ56Z12TmK4HtgN2B7TJzTmZ+o23RqZIdly5l2qpVAExbtYod\nly6tOSJJkiQNZ8wJeURsFRFXRMTBAJm5JjPvy8w17Q9PY7Xj0qU875xz2H3BAra68052X7CA551z\njkm5JEnSJDXmhDwznwQOr/JYTbzHjjySJ/bem63vuosXveUtbH3XXTyx9948duSRdYcmSZKkFqom\n1cuAg9oZiNrj6Z135v7zz99o2f3nn8/TO+9cU0SSJEkaSdWE/DTgfRFxQkS8ICKmRcQWjbd2BqnR\nm7ZqFbudeupGy3Y79dQNc8olSZI0uVRNnH8C7AVcCKwEngTWNdyebEt0GrPtrrlmwzSVFVdfvWH6\nynbXXFN3aJIkSWqhUtlD4Oy2RqG2+dUxxwDFXPKnd96Zez/zGba75poNyyVJkjS5jPnCQBGxFcXF\ngC7IzO9MSFQd4oWBpqYdly7d8IFj2qpVY/7A4YWBJElSp410YSCrrGjcOln33LKOkiSp21hlRePS\n6QTZso6SJKnbWGVF49LpBNmyjpIkqdtYZUXj0ukE2bKOkiSp21StsrIIGNvZoOpKwyXI937mMxOS\nlDeWdbz//PPZ7dRTN5R1tJKMJEmaiiol5Jm5sM1xaIrqdIJsWUdJktRtxlz2cCqIiIOA/1v+uhvw\ntcw8rcV6lj1sg/GWIew0yx5KkqROG6nsYeWEPCIOBM4EDgZ2BOZk5q0RcS7wncycFJeGjIhPA5dm\n5g0t2kzIe5AJuSRJ6rS21iEvO5wHfA94GfD5pn7WA8dV6bfdIuJZFB8UNknGJUmSpMmgapWVjwLf\nAPYDTm1quxWYNZ6ghkTE8RFxS0SsjYhLm9p2iogrIuLxiFgREUe36OJw4JvtiEWSJEmaCFUT8lnA\nxVnMd2me8/IwsMu4onrGfcA5wCUt2i4C1pbbOha4OCL2bVrn7cCX2xRLZZ28kqUkSZKmlqoJ+Vpg\n+jBtuwKPVux3I5l5ZWZeDaxuXB4R04GjgDMy8zeZuQy4CpjfsM6WwCsz88Z2xFKVl3qXJEnSSKom\n5DcCJ0fEtIZlQyPl7wOuG1dUm7cPsC4z725YdhvFFJohh3cgjs3yUu+SJEkaSdULA50JLKNIgi+j\nSMYXRMT5wCuAiS5hsS0w2LRsENhu6JeyystmK73Mnz+fmTNnMnPmTObMmcPcuXPbGujQlSxf9Ja3\nbFjmpd4lSZK6W39/P/39/QwMDDAwMDDiuuMpezgLOI+i7OE0iuoqNwCnZuaPKnU6/LbOAWZm5nvL\n318O3JiZ2zascxpwcGa+dQz9TnjZw2mrVrH7ggVsfdddG5Y9sffeE3YlS22eZQ8lSVKntb3sIUBm\n3pqZr6UYlX4BsH1mHtruZHwYdwJbRsReDcsOAG7vwLbHpPFKliuuvnrD9JXtrpkUZdolSZJUs6pT\nVjbIzLXA/W2IZRPlHPVnUYzAbxkRWwNPZeaaiLgcWBQR76eo+vJm4NUTEcd4eKl3SZIkjaTylJVO\niIizgLPYuLTi2Zm5KCJ2Ai4FjqAotfjBzPzSGPv3Sp09yCkrkiSp00aasjKpE/KJZkLem0zIJUlS\np03IHHJJkiRJ42dCLkmSJNVozAl5RGwVEVdExMETEZAkSZLUS8ackGfmkxRXwXR0XZIkSRqnqkn1\nMuCgdgYiSZIk9aKqdchPA66MiMeBK4EH2Lg0IZm5fpyxSZIkSV2v6gj5T4C9gAuBlcCTwLqG25Nt\niU6SJEnqclVHyBfRNCIuSZIkaewqJeSZubDNcUiSJEk9yUopkiRJUo0qJ+QRcWBEXB4RD0fEUxEx\nq1x+bkQc2b4QJUmSpO5VKSGPiHnA94CXAZ9v6mc9cNz4Q5MkSZK6X9UR8o8C3wD2A05tarsVmDWe\noCRJkqReUbXKyizgqMzMiGiutvIwsMv4wpIkSZJ6Q9UR8rXA9GHadgUerdivJEmS1FOqJuQ3AidH\nxLSGZUMj5e8DrhtXVJIkSVKPqDpl5UxgGXAbcBlFMr4gIs4HXgHMbk94kiRJUnerNEKembcBBwMP\nAqcDAZxQNh+SmT9rT3iSJElSd6s6Qk5m3gq8NiK2AWYAv8rMNW2LrEMWL17MnDlzmDt3bt2hSJIk\nqcv09/fT398/4jqR2VwkZfMi4h7gbeVIeXPb/sDVmfniMXfcYRGRy5cvrzsMddjg4CCzZzurSpIk\ndU5EkJnRqq3qSZ17AlsP07YNsEfFfiVJkqSeUjUhh2eqqjR7JfCrcfQrSZIk9YxRzyGPiFOAU8pf\nE/hqRDzZtNqzKeaTf7E94UmSJEndbSwndd4DfKu8vwD4AfDLpnWeAO4APjX+0CRJkqTuN+qEPDOv\nAq6CYlI6sCgzV0xQXJIkSVJPqFT2MDP/uN2BSJIkSb2o8kmdEXFgRFweEQ9HxFMRMatcfm5EHNm+\nECVJkqTuVSkhj4h5wPeAlwGfb+pnPXDc+EOTJEmSul/VEfKPAt8A9gNObWq7FZg1nqAkSZKkXlFp\nDjlFwn1UZmZENNcjfxjYZXxhSZIkSb2h6gj5WmD6MG27Ao9W7FeSJEnqKVUT8huBkyNiWsOyoZHy\n9wHXjSsqSZIkqUdUnbJyJrAMuA24jCIZXxAR5wOvAGa3JzxJkiSpu1UaIc/M24CDgQeB04EATiib\nD8nMn7UnPEmSJKm7VR0hJzNvBV4bEdsAM4BfZeaatkUmSZIk9YDKCTlAROwO7A5sU/6+oS0za5tH\nHhF7ALcA/1UuentmrqorHkmSJGk4lRLyiHgxsBSYM7So/Jnl/QSmtXhoJ/Vn5jtqjkGSJEkaUdUR\n8k8BLwROBpYDT7YtovaZFxHXAzdm5ul1ByNJkiS1UrXs4WzgzzNzcWZem5nXN9/aEVxEHB8Rt0TE\n2oi4tKltp4i4IiIej4gVEXF0Q/P9wF6ZeQiwS0S8rR3xSJIkSe1WNSH/OZ0ZFb8POAe4pEXbRRQX\nKNoFOBa4OCL2BcjMdZn5m3K9K4ADOhCrJEmSNGZVE/JzgQ9GxHPaGUyzzLwyM68GVjcuj4jpwFHA\nGZn5m8xcBlwFzC/bt21Y/TXAXRMZpyRJklRVpTnkmbkkIl4GDETETcAjm66SC8Yd3fD2AdZl5t0N\ny24DDinvz4uIjwC/BlYAZ0xgLJIkSVJlVausvAf4a+BpYBabTl/J8YW1WdsCg03LBoHtADLzGuCa\n0XQ0f/58Zs6cycyZM5kzZw5z585tb6SSJEnqOf39/fT39zMwMMDAwMCI61atsnI2xdzs92Xmryr2\nMR6PA9s3LdsBeGysHS1ZsqQtAUmSJElD+vr66Ovr2/B74/V6mlWdQ74zcFFNyTjAncCWEbFXw7ID\ngNtrikeSJEmqpGpCfiOwbzsDaSUipkXENhQXGdoyIraOiGmZuQa4HFgUEdMjYh7wZsDhbkmSJE0p\nVRPyk4BPm89oAAAgAElEQVT3R8QxEbFzRGzRfGtTfGcAa4APAseU94cu8nM8MB14CPgccFxm/rRN\n25UkSZI6IjLHfv5lRKwv7w734MzMqvPTOyYicvny5XWHoQ4bHBxk9uzZdYchSZJ6SESQmS0nkldN\nmhcx8ZVUJEmSpK5XtQ75wjbHIUmSJPWkds31liRJklTBuOZ5R8QBwEuBbZrbMvOz4+lbkiRJ6gVV\nr9S5I/B14KChReXPxnnlJuSSJEnSZlSdsnIuxcWBDqZIxt8GHAYsBe4B5rQlOkmSJKnLVU3IX0+R\nlN9U/v7zzOzPzHcD36SoUy5JkiRpM6om5LsCKzLzaWAtsF1D2+XAG8cbmCRJktQLqibkvwBmlPdX\nAq9qaNt7XBFJkiRJPaRqlZUbKU7ovApYApwVEXsCTwELgKvbEZwkSZLU7aom5GcDu5X3z6M4wfOd\nwHSKZPzE8YcmSZIkdb+qV+q8G7i7vL8OOK28SZIkSRqDMc8hj4itIuKKiDh4IgKSJEmSesmYE/LM\nfBI4vMpjJUmSJG2salK9jGeu0ilJkiSpoqondZ4GXBkRjwNXAg8A2bhCZq4fZ2ySJElS16s6Qv4T\nYC/gQoo65E8C6xpuT7YlOkmSJKnLVR0hX0TTiLgkSZKksata9nBhm+OozeLFi5kzZw5z586tOxRJ\nkiR1mf7+fvr7+0dcJzKrD3RHxAHAS4FtmpoyM5dU7rhDIiKXL19edxjqsMHBQWbPnl13GJIkqYdE\nBJkZrdoqjZBHxI7A14FXUUxdGeq8Mbuf9Am5JEmSVLeqJ3WeC+wMvIYiGX8bcBiwFLgHmNOW6CRJ\nkqQuVzUhfz1FUn5T+fvPM7M/M98NfBM4qR3BSZIkSd2uakK+K7AiM58G1gLbNbRdDrxxvIFJkiRJ\nvaBqQv4LYEZ5fyXFXPIhe48rIkmSJKmHVK1DfiNwEHAVxcmbZ0XEnsBTwALg6nYEJ0mSJHW7qgn5\n2cBu5f3zKE7wfCcwnSIZP3H8oUmSJEndr+qFge4G7i7vrwNOK2+SJEmSxqDqHHJJkiRJbVB1ygoA\nEbE9sD8wE7gP+ElmPtaOwCRJkqReUDkhj4gPU0xT2Zbi4kAJPB4R52XmR9oUnyRJktTVKiXkEXE2\ncCbwKeCLwIPA84CjgbMjYsvMXNiuICVJkqRuVXWE/P3AxzLzLxuW3Q5cFxGPAh8AFo4zNkmSJKnr\nVT2pcwfgG8O0XVO2S5IkSdqMqgn5zcDsYdpml+2SJEmSNqPqlJU/B66IiKeAL/PMHPJ3AO8F3hoR\nG5L9zFw/3kDHIiKeC1wBrKO4eugxmflgJ2OQJEmSRiMyc+wPihhKsFs9OJqWZ2aOq7ziWEVEZPnE\nImIBMDMzz22xXi5fvryToWkSGBwcZPbs4b7gkSRJar+IIDOjVVvVRHkRrZPxSSE3/pSxHcUJp5Ik\nSdKkUykh71RJw4g4HngP8DvA5zPzvQ1tOwGXAkcAvwQ+lJlfaGg/APgkxQmmr+tEvJIkSdJYVT2p\ns1PuA84BLmnRdhGwFtgFOBa4OCL2HWrMzNsy8yCKeukf6kCskiRJ0phN6oQ8M6/MzKuB1Y3LI2I6\ncBRwRmb+JjOXAVcB88v2ZzWsPgj8ukMhS5IkSWPS0ZMt22gfYF1m3t2w7DbgkPL+yyPi7ykqrKyl\nqPwiSZIkTTpTNSHflmLku9EgxQmcZOYtPJOcj2j+/PnMnDmTmTNnMmfOHObOndveSCVJktRz+vv7\n6e/vZ2BggIGBgRHXnaoJ+ePA9k3LdgAeG2tHS5YsaUtAkiRJ0pC+vj76+vo2/B7RsuIhMMnnkI/g\nTmDLiNirYdkBWN5QkiRJU8y4EvKIOLRdgQzT/7SI2AaYRpGAbx0R0zJzDXA5sCgipkfEPODNgMPd\nkiRJmlI2O2UlIl5DcfXNVhYA325rRBs7AziLZy5CdAxwNsWFiY6nqEP+EPAwcFxm/nQCY5EkSZLa\nbjRzyPcF3g/8V4u2V7Y3nI1l5tkUCXirtkeAt03k9iVJkqSJttmEPDP/OSK2ysyPN7dFxJ9OTFiS\nJElSbxjtHPJPtVqYmRe3MRZJkiSp54wqIc/MtRMdiCRJktSLpmrZQ0mSJKkrjCkhj4gZEfE7EbH1\nRAUkSZIk9ZKxjpCfDDwFvDEi3hkRezavEBELImKrNsQmSZIkdb2xJuTPAV4GXJmZXwJ2i4jXNa1z\nBfBn7QhOkiRJ6nZjTcj/CpgNrIiIxRQXDPrPiHj70AqZOdjG+CRJkqSuNpoLA22QmeuAD0XEJ4ET\nga8A2wA/jIj/BXwTmAG8qN2BSpIkSd2oUpWVzFyZmX8B7Aq8DrgcWAkcB9xDMddckiRJ0maMaYS8\nWWYm8P3yRnmS55zMfGrckUmSJEk9YFwJebPMHAAG2tmnJEmS1M28MJAkSZJUIxNySZIkqUYm5JIk\nSVKNTMglSZKkGpmQS5IkSTWqnJBHxIERcXlEPBwRT0XErHL5uRFxZPtClCRJkrpXpYQ8IuYB3wNe\nBny+qZ/1FBcImhIWL17MzTffXHcYkiRJ6kL9/f0sXLhwxHWiuLbP2ETEjcAq4H8B04AngVdm5q0R\ncRTwD5n5wjF33GERkcuXL687DHXY4OAgs2fPrjsMSZLUQyKCzIxWbVUvDDQLOCozMyKaM/qHgV0q\n9itJkiT1lKpzyNcC04dp2xV4tGK/kiRJUk+pmpDfCJwcEdMalg2NlL8PuG5cUUmSJEk9ouqUlTOB\nZcBtwGUUyfiCiDgfeAXgBF1JkiRpFCqNkGfmbcDBwIPA6UAAJ5TNh2Tmz9oTniRJktTdqo6Qk5m3\nAq+NiG2AGcCvMnNN2yKTJEmSekDlhHxIZq4F7m9DLJIkSVLPqZSQR8RIJ22up6iy8kPgksx8sMo2\nJEmSpF5QdYQ8gH0oShyuoJhL/jzgRcAD5e9vAE6JiEMy8442xCpJkiR1naplD8+nqEX+yszcKzNf\nnZl7UVRXWQucDbwE+CXwN22JVJIkSepCVRPyjwALyxM7N8jMH1Ik4x/JzJ8D51FUY5EkSZLUQtWE\nfB+K0e9WfgnsXd6/G3hOxW1IkiRJXa9qQj4AvH+Ytg+U7QC/BayquA1JkiSp61U9qXMR8LmI+E/g\nK8BDwHOBPwD2B95Vrnc4cPN4g5QkSZK6VaWEPDO/EBEPU8wX/xDwLGAd8APgdZn5zXLVU4Gn2xGo\nJEmS1I3Gc6XOa4FrI2ILiqkpD2fm+qZ11o4zvkoiYnvgWmBf4CDLLkqSJGmyqjqHfIPMXJ+ZDzUn\n4zX7NUUd9MvqDkSSJEkaSeURcoCIOAB4KbBNc1tmfnY8fY9HZj4NrIqIqCsGSZIkaTQqJeQRsSPw\ndeCgoUXlz2xYbdwJeUQcD7wH+B3g85n53oa2nYBLgSMoSi1+KDO/MN5tSpIkSZ1UdcrKucDOFBf9\nCeBtwGHAUuAeYE5booP7gHOAS1q0XURxVdBdgGOBiyNi3zZtV5IkSeqIqgn56ymS8pvK33+emf2Z\n+W7gm8BJ7QguM6/MzKuB1Y3LI2I6cBRwRmb+JjOXAVcB81t047QVSZIkTVpVE/JdgRXlXO21wHYN\nbZcDbxxvYJuxD7AuM+9uWHYbsN/QLxHxdYrpLP8cEe+e4HgkSZKkSqqe1PkLYEZ5fyXwKqC//H3v\nccY0GtsCg03LBmn4YJCZo/pQMH/+fGbOnMnMmTOZM2cOc+fObWOYkiRJ6kX9/f309/czMDDAwMDA\niOtWTchvpDih8ypgCXBWROwJPAUsAK6u2O9oPQ5s37RsB+CxsXa0ZMmStgQkSZIkDenr66Ovr2/D\n7yMV/6uakJ8N7FbeP4/iBM93AtMpkvETK/Y7WncCW0bEXg3TVg4Abp/g7UqSJEltNeY55BGxFfD3\nlCdLZua6zDwtM1+QmTMy812ZuaodwUXEtIjYBphGkYBvHRHTMnMNxVz1RRExPSLmAW+mGK2XJEmS\npowxJ+SZ+SRweJXHVnAGsAb4IHBMef/0su14ihH5h4DPAcdl5k87EJMkSZLUNlWnrCyjmEPe375Q\nNpWZZ1NMj2nV9ghF/XNJkiRpyqqakJ8GXBkRjwNXAg+w8VU6ycz144xNkiRJ6npVp538BNgLuJCi\n7OGTwLqG25NtiU6SJEnqclVHyBfRNCIuSZIkaewqJeSZubDNcUiSJEk9qROVUiRJkiQNo3JCHhEH\nRsTlEfFwRDwVEbPK5edGxJHtC1GSJEnqXpUS8vJCPN8DXgZ8vqmf9cBx4w9NkiRJ6n5VR8g/CnwD\n2A84tantVmDWeIKSJEmSekXVKiuzgKMyMyOiudrKw8Au4wtLkiRJ6g1VR8jXUly2vpVdgUcr9itJ\nkiT1lKoJ+Y3AyRExrWHZ0Ej5+4DrxhWVJEmS1COqTlk5E1gG3AZcRpGML4iI84FXALPbE54kSZLU\n3SqNkGfmbcDBwIPA6UAAJ5TNh2Tmz9oTniRJktTdqo6Qk5m3Aq+NiG2AGcCvMnNN2yKTJEmSekDl\nhHxIZq4F7m9DLJIkSVLPqXphoB9FxMkR8bx2ByRJkiT1kqpVVh4AzgPujYh/j4g/KqeuSJIkSRqD\nqid1vgGYCfwfiosAfR54MCIujYhD2xifJEmS1NWqjpCTmQ9l5j9k5iuB/YBPAIcB34yIle0KcKIt\nXryYm2++ue4wJEmS1IX6+/tZuHDhiOtEZo64wmiVU1b+APgosFtmTtvMQ2oXEbl8+fK6w1CHDQ4O\nMnu2pfIlSVLnRASZGa3aKo+QN3R+WER8mqIm+WeBnwMnjrdfSZIkqRdUKnsYEfsDxwLvAl4ADAAX\nAksy87/bFp0kSZLU5arWIf9P4FHgy8BnM/PG9oUkSZIk9Y6qCfk7gK9m5hPNDRFxCLAgM987rsgk\nSZKkHlC17OFljcl4ROwdEYsiYgXwbYqEXZIkSdJmVD6pMyJ2iIgPRMQy4GfA6cAjwJ8Cu7UpPkmS\nJKmrjSkhj4gtIuINEfEliqt1/hOwB0UNcoCTM/OTmTnY5jglSZKkrjTqhDwiPgbcB3wVeBNwBXAk\n8ELgw0DLuopSN1i9ejU33HADjzzySN2hSJKkLjOWEfJTgOcC/wa8MDOPycz/yMz1QHuuLiRNQhdc\nsJRZsz5BX99aDjzw41xwwdK6Q5IkSV1kLAn5JcBjwBuBn0XExyNizsSEJU0Oq1ev5sIL72HlyjNZ\nv/4IVq48kwsvvJvVq1dP+HYdkZckqTeMOiHPzPcDzweOAX4A/G/gexHxU+CDOEquLnT77bdz770H\nbbTs3ntfxR133DFh23REXpKk3jKmkzozc21mfiEzh+aO/zXwNPBXFHPIPxoRx0bENu0PVeq8/fff\nn913v2mjZbvv/j3222+/CdleXSPykiSpPpXLHmbmA5n5d5m5PzCHotLKS4DPUlRgkaa8nXbaiZNO\nejF77LGILba4lj32WMRJJ+3FTjvtNCHbq2NEHpwiI0lSnSon5I0y8weZeSJF/fE/APrb0a80GZxy\nyjHceusJXH/9s/nRj07klFOOmbBtdXpEHpwiI0lS3dqSkA/JzHWZeUVmvq2d/Up1mzFjBvPmzZuw\nkfEhnR6Rd4qMJEn1i8zuOxczIrYHrgX2BQ7KzJbf90dELl++vKOxqX6Dg4PMnj277jBGtHr1au64\n4w7222+/Cf0QcMMNN9DXt5b164/YsGyLLa7l+uufzbx58yZkm6tXr+b2229n//33n/APOJIkTRYR\nQWa2vG5PW0fIJ5FfA28ALqs7EKmKTo3Id3qKjNNjJEnaVFcm5Jn5dGauwquHSiPq5BQZp8dIktTa\npEvII+L4iLglItZGxKVNbTtFxBUR8XhErIiIo+uKU+oWnTppta4KMpIkTXaTLiEH7gPOobgyaLOL\ngLXALsCxwMURsS9ARJwSEddFxGkdi1TqEp2YIlNHBRmwpKMkafKbdAl5Zl6ZmVcDG32PHRHTgaOA\nMzLzN5m5DLgKmF8+7oLMPCwzP9bUpdNWpEmg0xVkwDnrkqSpYdJWWYmIc4CZmfne8veXAzdm5rYN\n65wKHJKZb23x+K8DBwArgU9m5mdbrGOVlR40FaqsdLNOVZBZvXo1s2Z9gpUrz9ywbI89FnHrrScw\nY8aMCd2uVWQkSc1GqrKyZaeDGYdtgcGmZYPAdq1Wzsw3jqbT+fPnM3PmTGbOnMmcOXOYO3fuOMOU\nNJKh6TETbaQ56xO1/QsuWMqFF97DvfcexO67f5yTTnrxhF5ISpI0efX399Pf38/AwAADAwMjrjuV\nEvLHge2blu0APDaeTpcsWTKeh0uapIo56x9n5cpnaqwXc9ZPnJDtNVaRAVi58gguvHARCxasntAR\neUnS5NTX10dfX9+G3yOGn0U96eaQj+BOYMuI2Kth2QHA7TXFI2kS6/Sc9V6pIuNJspLUfpMuIY+I\naRGxDTCNIgHfOiKmZeYa4HJgUURMj4h5wJsBh7gltdSpko5QXxWZTvIkWUmaGJMuIQfOANYAHwSO\nKe+fXrYdD0wHHgI+BxyXmT+tI0hJU0OnrnpaRxWZTvLCTpI0cSZtlZVOsMpKb7LKiiZSp6rIdNoN\nN9xAX99a1q9/Zk7+Fltcy/XXP7sjJ+l2ilVyJE2UkaqsTMYRckmasjo1Ij+kU3O6nZIjSRPHhFyS\npqhOJpB1Tcnp1AeOuqbkeJKsJDAhl6QpqY4EspMnyUJnP3DUUSXHEXlJQ0zIJWkKqqvMYqem5HT6\nA0enp+T0yoi83wC0j/uyu5mQS9IU1O1zujv9gaMX6tZ3ekTebwDax33Z/UzIJWkK6vYyi3V84Ojm\nuvWdHpHvhTKZnt+gdjIhl6QpqtNzujuprg8c3Vq3vtMj8t1+5VrPb1C7WYfcOuQ9xzrk0tTRrXXd\nh3Tq+T3yyCMceODHWbnyzA3L9thjET/60YkTst1Ob6+TVq9ezaxZn9jkud166wnMmDGj7dvr9L7s\n9PPrJdYhlyRNSZ2u695p3Toi381Tqjy/QRPBEXJHyHuOI+SSelWnv3Ho5PY6dZXVukb/u/XblF7i\nCLkkSer4Nw6d2l4vXCSrW79NGdLrJ5E6Qu4Iec9xhFySukddc549v6F9LrhgKRdeeA/33nsQu+9+\nEyed9OKuOkl9yEgj5CbkJuQ9x4RckrrHDTfcQF/fWtavP2LDsi22uJbrr3828+bNqzEyjUYvnUTq\nlBVJktSVuv0iWd3Ok0gLJuSSJGnK6uaKLr3AD1QFp6w4ZaXnOGVFkrpPt8/p7mbFHPK7uffeV7H7\n7t/jpJP2cg55LzEh700m5JIkTS698IFqpIR8y04HI0mSJDUaKuvYq3p+DvnixYu5+eab6w5DkiRJ\nHdLJuuf9/f0sXLhwxHWcsuKUlZ7jlBVJknpXXXXPLXsoSZKknrd69WouvPAeVq48k/Xrj2DlyjO5\n8MK7Wb169YRvdyQm5JIkSeoJddQ9v+CCpcya9YkR1zEhlyRJUk/odN3zxhH5kZiQS5IkqSd0+kJS\nrUbkW/GkTk/q7Dme1ClJUm/rVN3zRx55hAMP/Hg5Qu5JnZIkSRLwTN3zib4IUeOI/EgcIXeEvOc4\nQi5Jkjpp9erV7Lzzzo6QS5IkSXWYMWPGiO0m5JIkSVKNTMglSZKkGpmQS5IkSTUyIZckSZJqZEIu\nSZIk1ciEXJIkSarRlnUHMBEi4rnAFcA64CngmMx8sN6oJEmSpE116wj5LzPz9zKzD1gCvK/meCRJ\nkqSWujIhz40vP7odcHtdsUiSJEkjmXQJeUQcHxG3RMTaiLi0qW2niLgiIh6PiBURcfQI/RwQETcB\nxwO3TnTcUrv19/fXHYIkTSm+b2qqmnQJOXAfcA5wSYu2i4C1wC7AscDFEbEvQEScEhHXRcRpAJl5\nW2YeBJwJfKgjkUtt5D8WSRob3zc1VU26hDwzr8zMq4HVjcsjYjpwFHBGZv4mM5cBVwHzy8ddkJmH\nZebHIuJZDQ8dBH7dofBrc/PNN3fFNsfbZ9XHj+VNfDTrtmudblDX82z3dtvRX5U+2n1sjmY9j82p\nt93x9jnRx+Zo1/e98xl1PM9uOTbH+rjJ8N456RLyEewDrMvMuxuW3Qbs12Ldl0fE9RHxLeAk4LxO\nBFin73//+12xzfH2WfXxJuQTp1uSHhPy7tMtx2Y7+jQhn3xMyMf3+Kn23hkbn/84eUTEOcDMzHxv\n+fs84F8zc7eGdf4EeFdmHlZxG5PzyUuSJKnrZGa0Wj6V6pA/DmzftGwH4LGqHQ63UyRJkqROmUpT\nVu4EtoyIvRqWHYAlDSVJkjSFTbqEPCKmRcQ2wDSKBHzriJiWmWuAy4FFETG9nMLyZooL/0iSJElT\n0qRLyIEzgDXAB4Fjyvunl23HA9OBh4DPAcdl5k/rCFKSJElqh0l7UqckSZLUCybjCLmkFiLiuRGx\nLCL6I+KbEfG8umOSpMkuIvaIiIfKiwdeFxE71x2T1MwRcmmKiIjI8g82IhZQlAU9t+awJGlSi4g9\ngPMy8x11xyINxxFyaYrIjT89b4cVhiRptOaVFwz8m7oDkVoxIZc6LCKOj4hbImJtRFza1LZTRFwR\nEY9HxIqIOLqp/YCIuIniBOdbOxm3JNVpHO+d9wN7ZeYhwC4R8baOBi6Nggm51Hn3AecAl7RouwhY\nC+wCHAtcHBH7DjVm5m2ZeRBwJvChDsQqSZNFpffOzFyXmb8p17uC4hom0qRiQi51WGZemZlXA6sb\nl0fEdOAo4IzM/E1mLgOuAuaX7c9qWH0Q+HWHQpak2o3jvXPbhtVfA9zVoZClUduy7gAkbbAPsC4z\n725YdhtwSHn/5RHx98BTFCNB7+1wfJI0GW3uvXNeRHyEYhBjBcX1TqRJxYRcmjy2pRj5bjRIcQIn\nmXkLz/yDkSQVNvfeeQ1wTaeDksbCKSvS5PE4sH3Tsh2Ax2qIRZKmCt87NeWZkEuTx53AlhGxV8Oy\nA7C8oSSNxPdOTXkm5FKHRcS0iNgGmEbxT2TriJiWmWuAy4FFETE9IuYBbwaW1BmvJE0Gvneqm5mQ\nS513BrAG+CBwTHn/9LLteGA68BDwOeC4zPxpHUFK0iTje6e6Vmx88T9JkiRJneQIuSRJklQjE3JJ\nkiSpRibkkiRJUo1MyCVJkqQamZBLkiRJNTIhlyRJkmpkQi5JkiTVyIRckiRJqpEJuaSuFRELImJ9\nRLy47liaRcQBEXFWROzYom19RCya4O1Pj4glEfFgub3zW6zz6bJtpNvTEXFwRBxS/n7wRMY9WfTa\n85U0sbasOwBJmmCT9XLELwfOApYAv6ph+8cD7wT+GPhv4IEW6ywCLm74/f3Ae4HfA9Y3LL+j/HlQ\nw/1u90N66/lKmkAm5JJUj6DeDwu/DdyfmUuHWyEzVwArhn6PiN8v734/M9e3eMj32xvi5BMRWwCR\nmY/TA89XUmc4ZUVSzyunH3wzIgYj4vGIuCYi9mtapz8ibvj/7d1/rNdVHcfx50tSQW3NlAubCCm4\nWlF/6FaCrdqKKInMH21h+WMtWXNqW+YwWTKmMZMyVm41KVxqzjEjppUGVxBc/hiGglkieCeWRuDN\nEaGAwLs/3ufLPvdzv3zvvQled7+vx/bZl8/5nHvO+ZwPjPc99/05V9KnJf1Z0k5Jz0j6UpP2Zkj6\nm6Q3JK2TNL18/Ypy/RJgUam+qZL6MbbWzpWSusq4Hpb0wX7ez9ckPV363ybpDkmjK9f3AxcDY6tp\nJwObtV599krhqMzZVElPSXpd0lpJH5U0TNI8Sa9I6i7pMSNqbY6Q9IMyB7vL53WS1MdYhkm6QdKm\nyhysljS5Vm9mbZ5+Ien4Wp39km6UNEtSF7AbmHiwlBVJ50l6rPz9eE3SYkkn1+pcWOZhh6TtktZL\numxgM25mQ4kDcjNra5KmAZ3Af4CvAjOAdwOPSDqpUjWA8cAC4IfAuWSax+JqjrqkKcBdZCrDuaXu\nAuC0Slu/B24sfz6fTH2YRM+0kYuAs4GrgEuBscDSskLb6n5mAncAz5b+ZwFTgYclHVOqnQksK/19\nrPS9tlW7/VRf8Q9gAnAzMA+4ADgauI9MhRkFXALMJed+TuU+hpUxfh34MfA5YCHwvdJeK9cC3yLn\n/bPk/D0EvLfS/k3AraWP6cB3Sh9/aBLwX0o+i6uBacArze5X0jeBe4G/kM91JjCRnPtjS52Pk2lK\nK4FzSr3bgF7vEphZG4kIHz58+BiSBxns7QNObVFnI7CsVnYcsA24pVK2klwdPbVSNhLYC1xbKXsU\nWF9r73Qy53pFf8ZW6m4AhlXKzi/1z2xxL0cAW4DOWnkj5/uKStmdQNcA53NOGcMRTa59slz7RJM5\nG1cpm17GUp/z3wAvVM4vKu2dVat3HbALOLHFOO8H7m1xfVx5brNr5ZPK2L5Yexb/AI5qdb/AseS7\nAAub9LUbuKqcXw28Otj/Nnz48PHOOrxCbmZtS9IEctX77pLmMKyszO4CHgPqaRwbI6KrcRIR24Ct\n5Op1I7/4DDK4pFJvLZVc7H5aHhH7KufPkHnnYw9SH+D9QAdwd63/PwGbySDy7fZ8RGyunD9XPv9Y\nq/ccMKZyPpUc8+O1Z7McOIpc5T+YNcDZJdXkLElH1q5PIeey/tzXADvo/dwfjIg9fdznJPInK/U2\nXy731mhzDXC8coebaZLe00e7ZtYGHJCbWTvrKJ+/BN6sHHvI1IQTavX/3aSN3cDw8ucTgSPJIL3u\nXwMcW72v3eVzeL1iRSMlo9mOKVsq199Or9XO97Qof1clJacDeB89n8ubwBNkqkj92VR9n1zNnw6s\nBrolLZLUuP8OMiB/gd7P/bgmbTebz7pGmw81aXNio82IWA18mfzmYwmwTdJySR/uRx9mNkR5lxUz\na2fd5fO7ZB55XV+ronWvkkFYR5Nro8gV38OpEcSPbnJtNPDkYe7/UOoGusjgtdlLnC8e7AvLTxbm\nA/MldQBfIPPQR5DvCHSTQf0Umm852V07789uOI2vuZjmWyHuqIxvCbCk5PR/isyJf4CePyEwszbi\ngED6XYAAAAKdSURBVNzM2lZEbJD0IvChiOjrRcH+tLdf0pNkvvfcRrmkM4BT6BmQN1a8e+ws8hZt\nIFfivwLcXul/MpnLPP8Q9nW4PQicB+yMiOf/30YiYiuwqLy8O7EULydzw8dFxIq3PNL0KBl0nxYR\nd/VzbK+TL5GOBxZIOiEi6t8MmFkbcEBuZkOdgM9L2lIr3x4RneQvyFkq6WhgMbnKPQqYDGyOiAUD\n7G8OsEzSb8ndM0aWsn/S+5fpCLhC0q/IlfV1EbF3gP0dUL4huB74uaQ7yd1expA7umygEqQfJi23\nIxygX5O7m6yQ9CNgHZk7PoFMRTknInY1HYS0tNRfS6bGnE7uoPIzgIjoknQzcKukDwCryPcGxgKf\nIV/MXNWPMR6434jYIema0mYHueK9HTiJzN1fGRH3SJpL/v1aSe7WcjK5k85TDsbN2pcDcjMb6gL4\nSZPyZ4GPRMQDZS/p2eS2eiPIfOvHgXuatNWs/QPlEdEp6UIyCF8CbAK+Xc63V+qtlzSH3BrvG+Q7\nPacAL9Xb7KP/nhUiFkraCVwDLAX+S26zOCsi3hhoewMcQ3/HfLA2qvO4V9JUcgvDy8i52Unmff+O\n1ulEq8hUl8uBY8g5vYncerHR/mxJfyW/Ibu89P13Mgd8Y21MfY63tHmbpJfIuZ9B/h/7MvAI8HSp\n9gRwJXALmdO/lXzB9foW92NmQ5wi3qm/VdrMbGiQNIYM8m6IiHl91Tczs/bigNzM7BCSNJxc/ewk\n01/GkyumI4GJETHQ3VbMzGyIc8qKmdmhtY/c0eSn5FZ3O8mt9y5wMG5mZs14hdzMzMzMbBD5FwOZ\nmZmZmQ0iB+RmZmZmZoPIAbmZmZmZ2SByQG5mZmZmNogckJuZmZmZDSIH5GZmZmZmg+h/Pqx2QpxX\naOkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb30acaa7d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# pylab.rcParams['figure.figsize'] = (3.5,1.7)\n",
    "err2 = [np.mean(j) for _,j in err.items()]\n",
    "min_len = np.max(np.where(np.array(err2) >= 0.5))+1\n",
    "\n",
    "loglog(lengths[0:min_len], err2[0:min_len], 'x', color = 'r', mew=2, ms=5)\n",
    "loglog(lengths[min_len:], err2[min_len:], 'o', color = 'b', ms = 5)\n",
    "yticks([10**-3,10**-1,10**1,10**3,10**5, 10**7, 10**9], fontsize = 12)\n",
    "xticks([10**3,10**5], fontsize = 12)\n",
    "\n",
    "pareto_front = lengths[min_len]/10**(1.0/12)\n",
    "axvspan(100,pareto_front, alpha=0.3, color='gray')\n",
    "xlabel('Length of Time series', fontsize = 16)\n",
    "ylabel(r'Average $\\ell^1$ parameter error', fontsize = 16)\n",
    "title('Parameter error for diffusion equation', fontsize = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
